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Abstract. It is now generally assumed that the heterogeneity of most networks in 
nature probably arises via preferential attachment of some sort. However, the origin 
of various other topological features, such as degree-degree correlations and related 
characteristics, is often not clear and attributed to specific functional requirements. 
We show how it is possible to analyse a very general scenario in which nodes gain 
or lose edges according to any (e.g., nonlinear) functions of local and/or global 
degree information. Applying our method to two rather different examples of brain 
development - synaptic pruning in humans and the neural network of the worm C. 
Elegans - we find that simple biologically motivated assumptions lead to very good 
agreement with experimental data. In particular, many nontrivial topological features 
of the worm's brain arise naturally at a critical point. 
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1. Introduction 

The conceptual simplicity of a network - a set of nodes, some pairs of which connected 
by edges - often suffices to capture the essence of cooperation in complex systems. 
Ever since Barabasi and Albert presented their evolving network model [1], in which 
linear preferential attachment leads asymptotically to a scale-free degree distribution 
(the degree, k, of a node being its number of neighbouring nodes), there have been 
many variations or refinements to the original scenario [21 El IH El El [7] (for a review, 
see [8]). In [9], we showed how topological phase transitions and scale-free solutions 
could emerge in the case of nonlinear rewiring in fixed-size networks. Now we extend 
our scope to more general and realistic situations, considering the evolution of networks 
making only minimal assumptions about the attachment /detachment rules. In fact, 
all we assume is that these probabilities factorize into two parts: a local term which 
depends on node degree, and a global term, which is a function of the mean degree of 
the network. 

Our motivation can be found in the mechanisms behind many real-world networks, 
but we focus, for the sake of illustration, on the development of biological neural 
networks, where nodes represent neurons and edges play the part of synaptic interaction 
flU\ [TTl [T2] . Experimental neuroscience has shown that enhanced electric activity 
induces synaptic growth and dendritic arborization [131 HH UHl HE] . Since the activity of 
a neuron depends on the net current received from its neighbours, which is higher the 
more neighbours it has, we can see node degree as a proxy for this activity - accounting 
for the local term alluded to above. On the other hand, synaptic growth and death 
also depend on concentrations of various molecules, which can diffuse through large 
areas of tissue and therefore cannot in general be considered local. A feature of brain 
development in many animals is synaptic pruning - the large reduction in synaptic 
density undergone throughout infancy. Chechik et al. [TTl [18] have shown that via 
an elimination of less needed synapses, this can reduce the energy consumed by the 
brain (which in a human at rest can account for a quarter of total energy used) while 
maintaining near optimal memory performance. Going on this, we will take the mean 
degree of the network - or mean synaptic density - to reflect total energy consumption, 
hence the global terms in our attachment/detachment rules. 

An alternative approach would be to consider some kind of model neurons 
explicitely and couple the probabilities of synaptic growth and death to neuronal 
dynamic variables, such as local and global flelds. In a Hopfleld network, for example, 
the expected value of the fleld (total incoming current) at node i is proportional to 
i's degree [19], the total current (energy consumption) in the network therefore being 
proportional to the mean degree; qualitatively, these observations are likely to hold 
also in more realistic situations [20], although relations need not be linear. Co-evolving 
networks of this sort are currently attracting attention, with dynamics such as Prisoner's 
Dilemma [21], Voter Model [22] or Random Walkers [23]. Although we consider this line 
of work particularly interesting, for generality and analytical tractability we opt here 
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to use only topological information for the attachment/detachment rules, although our 
results can be applied to any situation in which the dynamical states of the elements at 
the nodes can be functionally related to degreeiji. 

Following a brief general analysis, we show how appropriate choices of functions 
induce the system to evolve towards heterogeneous (sometimes scale-free) networks while 
undergoing synaptic pruning in quantitative agreement with experiments. At the same 
time, degree-degree correlations emerge naturally, thus making the resulting networks 
disassortative - as tends to be the case for most biological networks - and leading to 
realistic small-world parameters. 

2. Basic considerations 

Consider a simple undirected network with nodes defined by the adjacency matrix 
a, the element aij representing the existence or otherwise of an edge between nodes i 
and j. Each node can be characterized by its degree, ki = Initially, the degrees 

follow some distribution p{k,t = 0) with mean K{t). We wish to study the evolution 
of networks in which nodes can gain or lose edges according to stochastic rules which 
only take into account local and global information on degrees. So as to implement 
this in the most general way, we will assume that at every time step, each node has a 
probability of gaining a new edge, p^^^^^ to a random node; and a probability of losing 

a randomly chosen edge, P^^. We assume these factorize as Pj = u{K)T:{ki) and 
P}'^^^ = d{K)a{ki), where u, d, vr and a can be arbitrary functions, but impose nothing 
else other than normalization. 

For each edge that is withdrawn from the network, two nodes decrease in degree: 
i, chosen according to cr{ki), and j, a random neighbour of i's; so there is an added 
effective probability of loss kj/{K,N). Similarly, for each edge placed in the network, 
not only / chosen according to 7r{ki) increases its degree; a random node m will also 
gain, with the consequent effective probability A^~^ (though se^. Let us introduce the 
notation 7t{k) = TT{k) + A^~^ and 5"(fc) = cr(A;) + k/{nN). Network evolution can now 
be seen as a one step process [21] with transition rates u{K)Ti{k) and d{K)a{k). The 
expected value for the increment in a given p(fc, t) at each time step (which we equate 
with a temporal derivative) defines a master equation for the degree distribution |9]: 



- [u (/t) 7r(A;) + d (k) cr{k)] p{k, t). 

Assuming now that p{k,t) evolves towards a stationary distribution, Pst{k), then 
this must necessarily satisfy detailed balance since it is a one step process [21]; i.e., the 

I For instance, the stationary distribution of walkers used for edge dynamics in [23] is actually obtained 
purely from topological information, although it can only be written in terms of local degrees for 
undirected networks. 

§ We are ignoring the small corrections that arise because j ^ i and I ^ to, which in any case would 
disappear if self-connections were allowed. 




(1) 
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flux of probability from A; to A; + 1 must equal the flux from k + 1 to k, for all k 
This condition (sufficient for ([1]) to be zero) can be written as 



Pstik), (2) 



dk L'^('*st) ^(^ + 1) 
where we have substituted a difference for a partial derivative and Kgt = Sfc^Pst(^)- 
Setting 71 and a so as to be normalized to one (i.e., '^kP{k)7i{k) = '^kP{k)a{k) = 1, Vt), 
which is equivalent to saying that at each time step exactly m(k) nodes are chosen to gain 
edges and d{K) to lose them, then in the stationary state we will have u{ngi) = d^Kg^) 
since the total number of edges will be conserved. From ([2]) we can see that Pst{k) will 
have an extremum at some value kg if it satisfies 7r(fce) = d'{ke + l). kg will be a maximum 
(minimum) if the numerator in ([2]) is smaller (greater) than the denominator for k > k^, 
and viceversa for k < ke- Assuming, for example, that there is one and only one such 
ke, then a maximum implies a relatively homogeneous distribution, while a minimum 
means Pst{k) will be split in two, and therefore highly heterogeneous. More intuitively, 
if for nodes with large enough k there is a higher probability of gaining edges than of 
losing them, the degrees of these nodes will grow indefinitely, leading to heterogeneity. 
If, on the other hand, highly connected nodes always lose more edges than they gain, 
we will obtain quite homogeneous networks. From this reasoning we can see that a 
particularly interesting case (which turns out to be critical) is that in which 7r(/c) and 
a{k) are such that 

n{k) = a{k) = v{k), VA;. (3) 

According to ([2]), condition ([3]) means that for large k, dpst{k)/dk — > 0, and Pstik) 
flattens out - as for example a power-law does. 

The standard Fokker-Planck approximation for the one step process defined by ([H) 

= \W ^ t^^'^)^^^) + n{K)m] p{k, t)} 
+^ {[diK)aik) - uiK)nik)] pik, t)} . (4) 
For transition rates which meet condition (j3]), (j4]) can be written as: 

^^ = l[di'^)+ui^)]^,[v{k)p{k,t)] 

+ [d{K)~u{K)]-^[v{k)p{k,t)]. (5) 

Ignoring boundary conditions, the stationary solution must satisfy, on the one hand, 
v{k)pgl{k) = Ak + B, so that the diffusion is stationary, and, on the other, u{Kgf^) = 
(i(Kgt), to cancel out the drift. For this situation to be reachable from any initial 
condition, u{k,) and d{K,) must be monotonous functions, decreasing and increasing 
respectively. 
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3. Synaptic pruning 

As a simple example, we will first consider global probabilities which have the linear 
forms: 

«K^)] = T^ 1--^ and d[K{t)] = --^, (6) 

where n is the expected value of the number of additions and deletions of edges per 
time step, and Kjnax is the maximum value the mean degree can have. This choice 
describes a situation in which the higher the density of synapses, the less likely new 
ones are to sprout and the more likely existing ones are to atrophy - a situation that 
might arise, for instance, in the presence of a finite quantity of nutrients. Again taking 
into account that tt and a are normalized to one, summing over Pf' — pP^^ we find 
that the increment in ^(t) is dn{t)/dt = 2{u[K{t)] - d[n{t)]} = 2{n/N) [1 - 2K(t)/Kmax] 
(independently of the local probabilities). Therefore, the mean degree will increase 
or decrease exponentially with time, from to l^max- Assuming that the initial 
condition is, say, /«(0) = ^max, and expressing the solution in terms of the mean synaptic 
density - i.e., p{t) = K{t)N / {2V), with V the total volume considered - we have 

PW=Pf(l + e-*/T^), (7) 

where we have defined p£ = pit — t- oo) and the time constant for pruning is rp = p^N/n. 
This equation was fitted in figured] to experimental data on layers 1 and 2 of the human 
auditory cortej|JJ] obtained during autopsies by Huttenlocher and Dabholkar [26] . 

It seems reasonable to assume that the initial overgrowth of synapses is due to the 
transient existence of some kind of growth factors. If we account for these by including 
a nonlinear, time-dependent term g{t) = aexp(— t/rg) in the probability of growth, i.e., 
u[K{t),t\ = {n/N)[l — K(t)/Kniax + g(t)], leaving d[K(t)] as before, we find that p(t) 
becomes 



pit) = Pi 



(8) 



where to is the time at which synapses begin to form (t = corresponds to the moment 
of conception) and Tg is the time constant related to growth. The inset in figure [T] shows 
the best fit to the auditory cortex data. Since the contour conditions and (for (I8]))f 
tf) are simply taken as the value of the last data point and the time of the first one, in 
each case, the time constants rp and rg are the only parameters needed for the fit. 



4. Phase transitions 



The drift-like evolution of the mean degree we have just illustrated with the example 
of synaptic pruning is independent of the local probabilities 7i{k) and a{k). The effect 
of these is rather in the diffusive behaviour which can lead, as mentioned, either to 

II Data points for three particular days (smaller symbols) are omitted from the fit, since we believe 
these must be from subjects with inherently lower synaptic density. 
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Conceptual age (days) 

Figure 1. Synaptic densities in layers 1 (red squares) and 2 (black circles) of the 
human auditory cortex against time from conception. Data from |26| . obtained by 
directly counting synapses in tissues from autopsies. Lines follow best fits to ([7|. 
where the parameters were: for layer 1, rp = 5041 days; and for layer 2, rp = 3898 
days (for we have used the last data pints: 30.7 and 40.8 synapses/z^m'^, for layers 
1 and 2 respectively). Data pertaining to the first year and to days 4700, 5000 7300, 
shown with smaller symbols, where omitted from the fit. Assuming the existence of 
transient growth factors, we can include the data points for the first year in the fit 
by using This is done in the inset (where time is displayed logarithmically). The 
best fits were: for layer 1, rg = 151.0 and rp — 5221; and for layer 2, rg = 191.1 and 
Tp — 4184, all in days (we have approximated to to the time of the first data points, 
192 days). 

homogeneous or to heterogeneous final states. A useful bounded order parameter to 
characterize these phases is therefore m = exp (— ct^/k^) , where = (/c^) — is 
the variance of the degree distribution ((■) = N^^^-{-) represents an average over 
nodes). We will use mgt = limj_j.oo m{t) to distinguish between the different phases, since 

= 1 for a regular network and m^f^ — > for one following a highly heterogeneous 
distribution. Although there are particular choices of probabilities which lead to 
(E]), these are not the only critical cases, since the transition from homogeneous to 
heterogeneous stationary states can come about also with functions which never meet 
condition ([3]). Rather, this is a classic topological phase transition, the nature of which 
depends on the choice of functions [271 1211 [29] . 

Evolution of the degree distribution is shown in figure [2] for critical and supercritical 
choices for the probabilities, as given by MC simulations (starting from regular random 
graphs) and contrasted with theory. (The subcritical regime is not shown since 
the stationary state has a distribution similar to the ones at t = 10^ MCS in the 
other regimes.) The disparity between the theory and the simulations for the final 
distributions is due to the build up of certain correlations not taken into account in our 
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Figure 2. Evolution of tlie degree distributions of networks beginning as regular 
random graphs with k(0) = 20 in the critical (top) and supercritical (bottom) regimes. 
Local probabilities are a{k) — k/{{k)N) in both cases, and 7r(fc) — 2a{k) — N~-^ and 
7r(fc) ~ k'^/'^ / {{k^/'^)N) for the critical and supercritical ones, respectively. Global 
probabilities as in ([G]), with n = 10 and Kmax — 20. Symbols in the main panels 
correspond to p{k,t) at different times as obtained from MC simulations. Lines result 
from numerical integration of ([!]). Insets show typical time series of k and m. Light 
blue lines are from MC simulations and red lines are theoretical, given by ([7]) and ^ , 
respectively. N = 1000. 



analysis. This is because the existence of some very highly connected nodes reduces 
the probability of there being very low degree nodes. In particular, if there are, say, x 
nodes connected to the rest of the network, then a natural cutoff, = x, emerges. 
Note that this occurs only when we restrict ourselves to simple networks, i.e., with only 
one edge allowed for each pair of nodes. This topological phase transition is shown in 
figure [3l where m^^ is plotted against parameter a for global probabilities as in (|6]) and 
local ones 7i{k) ~ k" and o"(fc) ~ k. This situation corresponds to one in which edges 
are eliminated randomly while nodes have a power-law probability of sprouting new 
ones (note that power-laws are good descriptions of a variety of monotonous response 
functions, yet only require one parameter). Although, to our knowledge, there are not 
yet enough empirical data to ascertain what degree distribution the structural topology 
of the human brain follows, it is worth noting that its functional topology, at the level 
of brain areas, has been found to be scale-free with an exponent very close to 2 |30] . 

In general, most other measures can be expected to undergo a transition along 
with its variance. For instance, highly heterogeneous networks (such as scale-free 
ones) exhibit the small-world property, characterized by a high clustering coefficient, 
C ^ {k)/N, and a low mean minimum path, I ~ In(A^) [3T]. A particularly interesting 
topological feature of a network is its synchronizability - i.e., given a set of oscillators 
placed at the nodes and coupled via the edges, how wide a range of coupling strengths 
will result in them all becoming synchronized. Barahona and Pecora showed analytically 
that, for linear oscillators, a network is more synchronizable the lower the relation 
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Figure 3. Phase transitions in m^^ for 7r(fc) ~ fc" and a{k) ~ k, and u{k) and 
d{K) as in N — 1000 (blue squares), 1500 (red triangles) and 2000 (black circles); 
k(0) = Kmax — 2n — N/50. Corresponding lines are from numerical integration of ([T]). 
The bottom left inset shows values of the highest eigenvalue of the Laplacian matrix 
(red squares) and of Q = (black circles), a measure of unsynchronizability; 

N — 1000. The top right inset shows transitions for the same parameters in the final 
values of Pearson's correlation coefficient r (see section [5]) , both for only one edge 
allowed per pair of nodes (red squares) and without this restriction (black circles) . 



Q = - where Xn and A2 are the highest and lowest non-zero eigenvalues of 

the Laplacian matrix (Aj^ = 6ijki — ciij), respectively [32]. The bottom left inset in 
figure [3] displays the values of Q and obtained for the different stationary states. 
There is a peak in Q at the critical point. It has been argued that this tendency of 
heterogeneous topologies to be particularly unsynchronizable poses a paradox given the 
wide prevalence of scale-free networks in nature, a problem that has been deftly got 
around by considering appropriate weighting schemes for the edges [331 El] (see alsq^ll, 
and |36j for a review). However, there is no generic reason why high synchronizability 
should always be desirable. In fact, it has recently been shown that heterogeneity can 
improve the dynamical performance of model neural networks precisely because the fixed 
points are easily destabilised [37] (as well as conferring robustness to thermal fluctuations 
and improving storage capacity [19]). This makes intuitive sense, since, presumably, one 
would not usually want all the neurons in one's brain to be doing exactly the same thing. 
Therefore, this point of maximum unsynchronizability at the phase transition may be a 
particularly advantageous one. 

On the whole, we find that three classes of network - homogeneous, scale-free (at the 
critical point) and ones composed of starlike structures, with a great many small-degree 

% Using pacemaker nodes, scale- free networks have also been shown to emerge via rules which maximize 
synchrony ^35j. 
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nodes connected to a few hubs - can emerge for any kind of attachment/detachment 
rules. It follows that a network subject to some sort of optimising mechanism, such as 
Natural Selection for the case of living systems, could thus evolve towards whichever 
topology best suits its requirements by tuning these microscopic actions. 

5. Correlations 

Most real networks have been found to exhibit degree-degree correlations, also known 
as mixing by degree |39]. They can thus be classified as assortative, when the degree 
of a typical node is positively correlated with that of its neighbours, or disassortative, 
when the correlation is negative. This property has important implications for network 
characteristics such as connectedness and robustness |10l ST]. A useful measure of 
this phenomenon is Pearson's correlation coefficient applied to the edges [39l HH [8]: 
r = {[kik'i] — [ki]"^) / {[kf] — [ki]"^), where ki and k'l are the degrees of each of the two nodes 
pertaining to edge /, and [■] = {{k)N)^^ represents an average over edges; |r| < 1. 

Writing = ^iji')^ ^ "^^^ be expressed in terms of averages over nodes: 



is the mean degree of the neighbours of node i, knn{k) is its average over all nodes such 
that ki = k. Whereas most social networks are assortative (r > 0) - due, probably, 
to mechanisms such as homophily [39] - almost all other networks, whether biological, 
technological or information-related, seem to be generically disassortative. The top 
right inset in figure |3] displays the stationary value of r obtained in the same networks 
as in the main panel and lower inset. It turns out that the heterogeneous regime is 
disassortative, the more so the larger a. (Note that a completely homogeneous network 
cannot have degree-degree correlations, since all degrees are the same.) It is known that 
the restriction of having at most one edge per pair of nodes induces disassortativity 

113]. However, in our case this is not the sole origin of the correlations, as can also 
be seen in the same inset of figure [3], where we have plotted r for networks in which we 
have lifted the restriction and allowed any number of edges per pair of nodes. In fact, 
when multiple edges are allowed, the correlations are slightly stronger. 

To understand how these correlations come about, consider a pair of nodes 
which, at a given moment, can either be occupied by an edge or unoccupied. We will 
call the expected times of permanence for occupied and unoccupied states r° and r^, 
respectively. After sufficient evolution time (so that occupancy becomes independent 
of the initial stat^ll), the expected value of the corresponding element of the adjacency 
matrix, E{aij) = iij, will be 



Note that this will always happen eventually since the process is ergodic. 



_ {k){ekuk))-{ey 
{k){k^)-{kY ' 

where knn{k) is the mean nearest-neighbour-degree function; i.e., if k. 




(9) 
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If (p- .) is the probability that will become occupied (unoccupied) given that it 



is unoccupied (occupied), then rf- ~ ^/v^j ^"^^ ^ij ~ yielding 




Taking into account the probability that each node has of gaining or losing an edge, 
we obtaiifl: = u{{k))N-^[7r{ki) + 7r{kj)] and p'j = d{{k))[a{ki)/ki + a{kj)/kj]. 
Then, assuming that the network is sparse enough that ^ pfj (since the number 
of edges is much smaller than the number of pairs), and particularising for power-law 
local probabilities nlk) ~ fc" and cr{k) ~ k^, the expected occupancy of the pair is 

. pfj _u{{k)) (k^) ( k^ + kf 



di{k)){k-)N \kt' + k^-')- 



Considering the stationary state, when u{{k)) = d{{k)), and for the case of random 

deletion of edges, /3 = 1 (so that the only nonlinearity is due to a), the previous 
expression reduces to 

{k) 



, , (fcf + fc"). (10) 

(Note that this matrix is not consistent term by term, since 'Y^- iij ^ ki, although it is 
globally consinstent: iij = {k)N.) The nearest-neighbour-degree function is now 



ki^ ' ' 2{k^) 



3 



(a decreasing function for any a), with the result that Pearson's coefficient becomes 
1 - (A;2)2(r) 



More generally, one can understand the emergence of these correlations in the 
following way. For the network to become heterogeneous, we must have Ti{k) + A^~^ > 
<y{k) + k/{{k)N) for large enough k, so that highly connected nodes do not lose more 
edges than they can acquire (see section [2]). This implies that TT{k) must be increasing 
and approximately linear or superlinear. The expected value of the degree of a node i, 
chosen according to 7r(/cj), is then E{k.j) = N~'^ Tlk'^i^)^ ~ while that of its 

new, randomly chosen neighbour, j, is only E{kj) = {k). This induces disassortative 
correlations which can never be compensated by the breaking of edges between nodes 
whose expected degree values are N~^^j^a{k)k and {k'^) / {k) if a{k) is an increasing 
function. It thus ensues that a scenario such as the one analysed in this paper will 
never lead to assortative networks except for some cases in which cr{k) is a decreasing 
function - meaning that less connected nodes should be more likely to lose edges. 

* Again, we are ignoring corrections due to the fact that i is necessarily different from j. 
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Assortativity could, however, arise if there were some bias also on the node chosen to be 
z's neighbour, e.g. on the postsynaptic neuron - which is precisely what happens in most 
social networks, where individuals do not generally choose their friends, partners, etc. 
randomly. Although there seem to be other reasons for the ubiquity of disassortative 
networks in nature [3^, it is possible that the generality of the scenario studied here 
may also play a part. 

We can use the expected value matrix e to estimate other magnitudes. For example, 
the clustering coefficient, as defined by Watts and Strogatz [31], is an average over nodes 
of Ci, with Ci the proportion of i's neighbours which are connected to each other; so 
its expected value is E{Ci) = iji conditioned to j and I being neighbours of z's. This 
means that, on average, we can make the approximation that kj = ki = (fc„ri) = 
+ (A;°+i)(A;-^)]/(2(A;")). Substituting this value in and taking into 
account that one edge of j's and one of /'s are taken up by i, we have 

For a rough estimate of the mean minimum path (the minimum path between two nodes 
being the smallest number of edges one has to follow to get from one to the other), we 
can procede as in [35]. For a given node, let us define the number of nearest neighbours, 
Zi, next-nearest neighbours, Z2, and in general mth neighbours, Zm- Using the relation 
Zm = zi {z2/zi)"^^^ , and assuming that the network is connected and can be obtained 
in / steps, this yields 
I 

]zm = N. (13) 
1 

On average, Zi = {k) and Z2 = {k)[{l — C){knn) — 1] (since for each second nearest 
neighbour, one edge goes to the reference node and a proportion C to mutual 
neighbours). Now, if N zi and Z2 ^ zi, (|T3l) leads to 

+ '■'(^/(^■» (14) 



1 + E 



6. The C. Elegans neural network 

There exists a biological neural network which has been entirely mapped (although 
not, to the best of our knowledge, at different stages of development) - that of the 
much-investigated worm C. Elegans [361 E]. With a view to testing whether such a 
network could arise via simple stochastic rules of the kind we are here considering, we 
ran simulations for the same number of nodes, = 307, and (stationary) mean degree, 
(k) = 14.0 (in the simple, undirected representation of the network). Using the global 
probabilities given by (E]) and local ones n{k) ~ fc" and a{k) ~ k (as in figure [3]), we 
obtain a surprising result. Precisely at the critical point, a = — 1.35, there are some 
remarkable similarities between the biological network and the ones produced by the 
model. 
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Figure S] displays the degree distributions, both for the empirical network and for 
the average (stationary) simulated network corresponding to the critical point, while 
the top inset shows the mean-nearest-neighbour degree function knnik) for the same 
networks. Both p{k) and knn{k) of the simulated networks can be seen to be very 
similar to those measured in the biological one. Furthermore, as is displayed in table 
[H the clustering coefficient obtained in simulation is almost the same as the empirical 
one. The mean minimum path is similar though slightly smaller in simulation, probably 
due to the worm's brain having modules related to functions |17]. Finally, Pearson's 
coefficient is also in fairly good agreement, although the simulated networks are actually 
a bit more disassortative. It should, however, be stressed that the simulation results are 
for averages over 100 runs, while the biological system is equivalent to a single run; given 
the small number of neurons, statistical fluctuations can be fairly large, so one should 
refrain from attributing too much importance to the precise values obtained - at least 
until we can average over 100 worms. Table [1] also shows the values of C, I and r both 
as estimated form the theory laid out in section [5l and for the equivalent network in the 
configuration model [39j - generally taken as the null model for heterogeneous networks, 
where the probability of an edge existing between nodes i and j is kikj/{{k)N). It is 
clear that whereas the conflgurat ion-mo del predictions deviate substantially from the 
magnitudes measured in the C. Elegans neural network, the growth process we are here 
considering accounts for them quite well. It is interesting that it should be at the 

Table 1. Values of small- world parameters C and /, and Pearson's correlation 
coefficient r, as measured in the neural network of the worm C. Elegans |46| , and 
as obtained from simulations in the stationary state {t = 10^ steps) for an equivalent 
network at the critical point when edges are removed randomly - i.e., for a = 1.35 
and (3 = 1. N = 307, k^^ = 14.0; averages over 100 runs and global probabilities 
as in ([6]). Theoretical estimates correspond to (|12p . (IT4l) and (fTTj) applied to the 
networks generated by the same simulations. The last column lists the respective 
configuration model values: C and I are obtained theoretically as in |39| . while r, from 
MC simulations as in [43], is the value expected due to the absence of multiple edges. 
(See also figure ID) 





Experiment 


Simulation 


Theory 


Conflg. 


c 


0.28 


0.28 


0.23 


0.15 


I 


2.46 


2.19 


1.86 


1.96 


r 


-0.163 


-0.207 


-0.305 


-0.101 



critical point that a structural topology so similar to the empirical one emerges, since it 
seems that the brain's functional topology may also be related to a critical point [l9l HH] . 
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Figure 4. Degree distribution (binned) of the C. Elegans neural network (circles) ^5] 
and that obtained with MC simulations (line) in the stationary state {t = 10^ steps) 
for an equivalent network in which edges are removed randomly (/3 = 1) at the critical 
point (a = 1.35). N = 307, Kg^ = 14.0, averages over 100 runs. Global probabilities 
as in ([6|). The slope is for fc~^/^. Top right inset: mean-neighbour-degree function 
knn{k) as measured in the same empirical network (circles) and as given by the same 
simulations (line) as in the main panel. The slope is for Bottom left inset: rrig^ 

of equivalent network for a range of a, both from simulations (circles) and as obtained 
with (p. (See also table [TJ) 



7. Discussion 

With this work we have attempted, on the one hand, to extend our understanding 
of evolving networks so that any choice of transition probabihties dependent on local 
and/or global degrees can be treated analytically, thereby obtaining some model- 
independent results; and on the other, to illustrate how such a framework can be applied 
to realistic biological scenarios. For the latter, we have used two examples relating to 
rather different nervous systems: 

i) synaptic pruning in humans, for which the use of nonlinear global probabilities 
reproduces the initial increase and subsequent depletion in synaptic density in good 
accord with experiments - to the extent that nonmonotonic data points spanning a 
lifetime can be very well fitted with only two parameters; and 

ii) the structure of the C. Elegans neural network, for which it turns out that by only 
considering the numbers of nodes and edges, and imposing random deletion of edges and 
power-law probability of growth, the critical point leads to networks exhibiting many of 
the worm's nontrivial features - such as the degree distribution, small-world parameters, 
and even level of disassortativity. 

These examples indicate that it is not farfetched to contemplate how many 
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structural features of the brain or other networks - and not just the degree 
distributions - could arise by simple stochastic rules like the ones considered; although, 
undoubtedly, other ingredients such as natural modularity |17], a metric [50] or 
functional requirements [12] can also be expected to play a role in many instances. 
We hope, therefore, that the framework laid out here - in which for simplicity we have 
assumed the network to be undirected and to have a fixed size, although generalizations 
are straightforward - may prove useful for interpreting data from a variety of fields. It 
would be particularly interesting to try to locate and quantify the biological mechanisms 
assumed to be behind this kind of network dynamics. 
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